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ABSTRACT 

Planets embedded in optically thick passive accretion disks are expected to produce perturbations in 
the density and temperature structure of the disk. We calculate the magnitudes of these perturbations 
for a range of planet masses and distances. The model predicts the formation of a shadow at the 
position of the planet paired with a brightening just beyond the shadow. We improve on previous 
work on the subject by self-consistently calculating the temperature and density structures under the 
assumption of hydrostatic equilibrium and taking the full three-dimensional shape of the disk into 
account rather than assuming a plane-parallel disk. While the excursion in temperatures is less than 
in previous models, the spatial size of the perturbation is larger. We demonstrate that a self-consistent 
calculation of the density and temperature structure of the disk has a large effect on the disk model. 
In addition, the temperature structure in the disk is highly sensitive to the angle of incidence of stellar 
irradition at the surface, so accurately calculating the shape of the disk surface is crucial for modeling 
the thermal structure of the disk. 

Subject headings: planetary systems: formation — planetary systems: protoplanetary disks — radia- 
tive transfer 



1. INTRODUCTION 

Giant planets forming by core accretion need to have 
cores of 10-20 to be massive enoug h to accrete 
a gaseous envelope (iHubickvj et al.l l2005l ). This pre- 
dicts that sizeable planet embryos form before cir- 
cumstellar gas disks dissipate. These disks are typ- 
ically modeled as passive accretion disks, optically 
thick and gas-dominated. Their temperature struc- 
ture is st rongly dependent on heat i ng from stellar ir- 
radiation (|Chiang fc Goldreichlll997l : ICalvet et"aDll99ll : 
iD^Alessio et al.l ll998L fl999al 120011 ). In particular, the 
disk temperatures depend strongly on the angle of inci- 
dence of the stellar irradiation on the disk surface. 

While a substantial amount of work on numerical hy- 
drodynamic simulations of planets embedded in disks 
has been carried out, calculating the radiative transfer 
of stellar irradiation in conjunction with th ese simu la- 
tions is prohibitively difficult. Simulations by lBate et al.l 
([2003.) illustrate the effect of hydrodynamics on the disk 
structure around a small embedded planet, but use an 
isothermal equation of state. The effects of MH D tur- 
bulence Jiaye been studied by iPapaloizou et al.l f2004) 
and Oishi et ajj (|2007i ). but ag!: ain, assuming a simple and 
unrealistic equation of sta te. iKlahr fc Kiev! ([2006) and 
iPaardekooper fc Mellemal (|2006l . 120081 ) have made an ef- 
fort to include radiative transfer as well as hydrodynam- 
ics in their calculations, but do not include the effects of 
stellar irradiaton on their models. 

The work presented in this paper does not include hy- 
drodynamics, but rather focuses on the effects of stellar 
irradiation, which is a particularly important source of 
heating in the photospheres of disks. Since it is the pho- 
tospheres of optically thick disks that are observed, the 
effects of stellar irradiation at the surface of disks are 
an important consideration in predicting observations 
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of protoplanetary disks. A limitation of the hydrody- 
namic simulations described above is that in order to 
adequately model the densest regions of the gas they are 
necessarily limited in spatial range above the midplane. 
The photosphere and surface of a disk are often several 
scale heights above the midplane. The models presented 
here are complementary to detailed hydrodynamic sim- 
ulations for this reason. 

We are particularly interested in determining if the 
growing cores of giant planets produce effects that are 
observable. If these effects are observed, they would af- 
firm the core accretion model as the paradigm for gi- 
ant planet formation. While a fully-formed Jupiter-mass 
planet would produce a larger feature in a disk, it would 
reveal little about how it formed. 

Previous work on small planets embedded in optically 
thick gas disks indicates that sub- Jovian mass planet 
cores can perturb the disk enough to alter the temper- 
ature structure of the disk in the immediate vicinity of 
the planet, with co nsequences for further evolution of 
th e planet (iJang-Condell fc Sasselov 2003 . 2004 . 2005h . 
In lJang-Condeh fc SasseloT(|200a l20o¥ 7henceforth Pa- 
per I and Paper II, respectively), the planet is predicted 
to gravitationally compress the disk in the vertical direc- 
tion, creating a shadow paired with a bright spot, leading 
to temperature variations. However, those models were 
limited by being plane-parallel, despite the sizes of the 
simulation boxes used. In addition, while the density 
perturbations were calculated under hydrostatic equilib- 
rium, they were not calculated self-consistently with the 
temperature perturbations. 

In this paper, we improve upon the model presented 
in Papers I and II by iteratively calculating the density 
and temperature structure of the disk for self-consistency 
and eliminating the assumption of a locally plane paral- 
lel model. The paper is organized as follows: in ^we 
describe the basic model and the improvements we have 
made, in ^we elaborate in detail on the method we use 
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for modeling radiative transfer, in 21 we describe our re- 
sults for a different planet masses and distances, in ^ 
we discuss our results in comparison to previous models, 
and in ^we present our conclusions. 

2. MODEL CALCULATION 

We adopt the f ormal ism developed in 

iJang-Condell fc Sasselovl (|2QQa I2QQ4D (henceforth. 
Paper I & II, respectively) for calculating the effects on 
temperature of shadowing and brightening of stellar illu- 
mination at the disk's surface. These methods are based 
on those of I Gal vet et al.l (pj991) and ID'Alessio et al.l 
(jl998l . Il999af ) in a 1-D plane-parallel disk. 

The improvements made on the previous model include 
iteratively calculating the temperature and density struc- 
ture of the disk for self-consistency. The new model is 
not plane-parallel as the previous models were. It ac- 
counts for the variation of disk structure with radius and 
the curvature of the disk in azimuth. 

Each iteration proceeds as follows. Starting with the 
density structure of the disk, we calculate the temper- 
ature of the disk in radiative equilibrium with viscous 
heating and stellar irradiation as heating sources. We 
adapt methods developed in Papers I and II to calcu- 
late radiative transfer on the surface of a perturbed disk. 
Then, given the new temperature structure, we recalcu- 
late the density structure assuming vertical hydrostatic 
equlibrium. 

Each of these steps are described in detail below. The 
initial conditions, which assume azimuthal symmetry, are 
calculated iteratively until they reach stability. The com- 
putationally intensive nature of the calculations prohibits 
large numbers of iterations when a planet is added, so rel- 
atively few iterations are carried out for the planet in the 
disk. 

2.1. Model Parameters 

We use rotating cyclindrical coordinates (r, ^, z) 
throughout, with the star at the origin, the planet at 
(a, 0,0), and the z-axis aligned with the orbital angular 
momentum vector. We adopt parameters for the stellar 
mass, radius, and effective temperature of = IM©, 
= 2.6 i?0, and = 4280 K, corresponding to an 
age of 1 Myr (|Siess et al.ll200Q) . The mass accretion rate 
is M = lO~^M0yr~^, and the viscosity parameter is 
ay = 0.01, which are typical of T-Tauri type stars. 

The simulation box is centered on the planet at 
(a, 0,0). The size of the box is scaled relative to the 
Hill radius, rHiii = a(mp/3M*)^/^. The box spans from 
Tmin = a - 12.5rHiii to Tmax = a + 12.5rHiii in r, and 
from ^min = -12.5rHiii/a to ^max = 12.5rHiii/a in cj). 
The height of the box is set to Zmax ~ 2zs where Zs 
is the height of the unperturbed disk surface at r = a. 
We assume symmtery across the midplane, so Zmin = 0. 
The box is decomposed onto a grid of 100 x 100 x 100 
points equally spaced in r, ^, and 2:, with the grid slice 
at = 0min used as a placeholder for structure of the 
unperturbed disk. 

We use the opacities from lD'Alessio etahl (|200lf ) using 
a dust model with parameters amax = 1 nim, T = 300 K, 
and p = 3.5, assuming that the dust opacities are con- 
stant throughout the disk. The values for the opacities 
(in cm^g~^) are as follows: the Rosseland mean opac- 
ity is xr = 1.91, the Planck mean opacity integrated 



over the disk spectrum (300 K) is /^p = 0.992, and the 
Planck mean opacities integrated over the stellar spec- 
trum (4000 K) are hip = 1.31 for absorption alone and 
Xp = 5.86 for absorption plus scattering. The absorp- 
tion fraction is then aabs = i^p/Xp^ while the scattered 
fraction is a = 1 — o^abs- The Rosseland mean opacity is 
used to calculate the photosphere of the disk, and Xp is 
used to calculate the surface of the disk. 

2.2. Heating Sources 

The basic calculation for radiative h eating is described 
in detail in iJang-Gondell fc Sasselovl (|2003). Here, we 
summarize those methods. The two main heating sources 
are from viscous heating (F^) and stellar irradiation (F^), 
offset by radiative cooling. 



A = PXrctbT^. 



(1) 



where p is density, ap is the Stefan-Boltzmann constant 
and T is the local temperature of the gas. Under equi- 
librium conditions, 



F^ = A. 



(2) 



If Ty and give the equilibrium temperature for solely 
viscous heating or stellar irradiation, respectively, then 
it follows that the equilibrium temperature given both 
sources of heating is 



i4\l/4 



(3) 



2.2.1. Viscous Heating 

We assume that viscous flux is generated at the mid- 
plane and transported radiatively in a grey atmosphere 
so that 

nl/4 



T„ = 



(Td + 2/3) 



(4) 



where Td is the optical depth perpendicular to the disk 
using the Rosseland mean opacity, xr' '^d = XrP^z' . 
The viscous flux emitted at the photosphere F^,ph at a 
distance r for a star of mass and radius accreting 
at a rate Ma is 



SGM^Ma 
47rr^ 



1/2- 



(5) 



('Pringlel [l98ll ). Viscous heating above the surface is as- 
sumed to be negligible. 

2.2.2. Stellar Irradiation 

The amount of heating from stellar irradiation is cal- 
culated using methods developed in Papers I and II 
for a plane-parallel disk and adapted for a fully three- 
dimensional system. The full details of this calculated 
are rather lengthy, so this has been set out into ^ 

2.2.3. Differential Rotation 

The gas moves in streamlines past the planet at ap- 
proximately the Kepler ian rate, resulting in the shearing 
out of the hot and cold spots as material moves into 
and out of shadows and brightenings. The calculation 
is done in steady state, meaning that the gas is treated 
as a steady flow so that temperatures as a function of 
spatial position is constant even though the gas itself 
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is in motion. Given the bulk velocity along a particular 
streamline, the movement of the gas from one grid cell to 
the next is equivalent to a time step equal to the length 
of the grid cell divided by the velocity. 

The equilibrium temperature calculated in equation 
([2j) represents the total amount of heating due to disk 
viscosity and stellar irradiation. We assume that the spe- 
cific heat per unit surface area of the disk is /cS/m, where 
k is the Boltzmann constant and fh is the mean molec- 
ular weight of the gas, which we assume to be primarily 
molecular hydrogen. Then we can approximate the rate 
at which a parcel of gas radiatively heats or cools as 



dT 
m dt 



(6) 



If Teq < T, the parcel of gas cools, if Teq > T the parcel 
of gas heats. 

2.3. Density Profile 

To calculate the density profile for the disk, we assume 
hydrostatic equilibrium. 



1 dP 

p dz 



-9z 



(7) 



where p is the density, P is the pressure, and Qz is the 
z-component of gravity. We assume the ideal gas law, 
P = pkT/fh. In the absence of the planet, the sole con- 
tribution to the gravity is the star: 



GM^z 



(r2 



2)3/2- 



(8) 



We use this equation to calculate the initial conditions. 
When a planet is added, the gravity has an additional 
contribution from the planet, so equation ([8]) becomes 



GM^z 



GMpZ 



(r^ + ^2)3/2 (^rp2 _^ rp2 _ 2r Tp COS (f) + z^Y"!^ 



■ (9) 



Given a vertical temperature profile, we calculate the 
density profile by integrating equation ([7j) from the top 
of the simulation box down the midplane. For the per- 
turbed disk, we require conservation of the total surface 
density. In a standard viscous accretion disk, the surface 
density of the disk and temperature of the disk are cou- 
pled. Thus, for the initial disk we normalize the total 
integrated surface density 



E = 2 / pdz' 



(10) 



with the surface density given by a steadily accreting 
viscous disk 



M 



r 



1/2' 



(11) 



where Uy is the viscosity of the disk ([Pringld 
Il98lf). We adopt a s tandard Shakura-Sunyaev viscos- 
itv (|Shakura fc Sunva ev 1973) with = ayC^H where 
ay is a dimensionless parameter, cq is the sound speed at 
the midplane, and H is the thermal scale height, given by 
H = cq/Qk where Qk is the Keplerian angular velocity. 



3. 3D RADIATIVE TRANSFER OF STELLAR IRRADIATION 

In this section, we describe the details of radiative 
transfer of stellar irradiation on the surface of a disk 
perturb ed by a plane t . The se me thods are based on the 
work of ICalvet etld] (|l99lf ) and iD'Alessio et all (|l998f ) 
which calculated the vertical temperature structure of a 
locally plane-parallel disk model without perturbations. 
In Papers I and II, we developed a method for using their 
solutions to calculate radiative transfer on a perturbed 
locally plane-parallel disk. Here, we extend the formal- 
ism to a fully three-dimensional disk. 

For a plane-parallel medium, the optical depth at disk- 
temperature frequencies perpendicular to the surface (rd) 
is related to the line-of-sight optical depth to stellar ra- 
diation (ts) as Ts = Xp'^d/iXRf^) where p is the angle 
of incidence of the stellar radiation at the surface. As 
shown in Paper II, for a locally plane-parallel disk is 
given by 



B{rs) ki + C2e + cse 

TT 47r 

(12) 

where p is the cosine of the angle of inciden ce of s tellar 
irradiation at the surface of the disk, P = >/3aabs, and 
the stellar fiux incident at the surface is 

where {rs^(j)s^ ^s) are the coordinates at the surface. The 
remaining coefficients are 

_ & + 9,iXr/x*p 6(1 - Xfl/X^) (3 - /?^) 



/?2 



/32(3 + 2/?)(1+/?m) 



02^ 



C3- 



X*P 


^fJ-XR 


) 


-3m') 


pnp 


Xp . 


1 (1- 




PX*P 


^xr\ 


(2- 


f 3m)(3-/?') 


Kp 


XpP) 


/3(3 4 


-2/3)(l-/?V) 



(15) 



(16) 



We adapt these equations to a three-dimensional disk 
by dividing the surface into grid elements by and 
and numerically integrating the contributions from each 
surface element. 

3.1. Flux Contributions 

For a surface element SA located at the contribution 
to the ffux at point P = (r, 0, z) in the disk is propor- 
tional to B{t^^) as given in equation [T2t with the optical 
depth calculated as 

rf' = 2/Z+- j\*ppdl (17) 

where the integral is carried out along the line segment 
PS and V is the cosine of the angle between PS and the 
surface element. The 2/3 term on the left hand side of the 
equation accounts for the optical depth from the star to 
the surface. The second term on the left hand side of the 
equation is what Tg would be between P and the infinite 
plane represented by the surface element. Defining dO. 
to be the solid angle subtended by the surface element, 
then the total heating from stellar irradiation at P in the 
limit as dO. ^ is 



Stot(P) 



(18) 
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For a plane parallel disk, this equation reverts to equa- 
tion [121 To calculate this numerically over an assemblage 
of finite surface elements, we assume that B(t^^^ii) 
stays fairly constant over a given surface element, and 
equation [TH] becomes 



(19) 



The expression for v d£t integrated over an arbitrary 
triangle with respect to P can be determined analytically 
as follows. We define si, S2, and S3, to be the vectors that 
point from P to the vertices of the triangle. Let n be the 
unit normal to the triangle and Oij be the angle begin 



the vectors and Sj. Then 



712 



smfyi2 



■(Si X S2) 



^23 



sm e/23 



■(S2 X S3) 



+ ^-^ S3 X Si 

sin6>3i 



(20) 




Fig. 1. — Angles used in calculating optical depths 



3.2. Surface Decomposition 

We define the surface of the disk to be where the optical 
depth to stellar irradiation integrated along the line of 
sight is Ts = 2/3. The optical depth at a given point 
P = (r, ^, z) in the disk is 



Xppdl, 



(21) 



integrating along the line segment which extends from 
the point toward the star, ending at the inner boundary 
of the box, at (rmin, ^in) where z^^ = zr^i^/r. The 
optical depth at the inner edge of the simulation box. 



Tin, is assumed to be 



(22) 



The density increases monotonically toward the mid- 
plane, so for every coordinate pair (ri,0j) the height of 
the surface hij can be uniquely determined. The surface 
is defined by the set of points Si^j = {vi^ hij). 

The set of neighboring points, Sij, 
and s^^j+i, define a surface element over which radia- 
tive transfer from stellar irradiation is calculated. The 
midpoint of this set of four points is 



cos 

h 



2 



^+i,j+i 



For the unperturbed disk hij is independent of j, and 
Sij^ Si-^ij, Si+i^j+i, and s^^j+i are co-planar. If n is the 
normal to the plane through these four points, then the 
cosine of the angle of incidence is 



M = Sri 



37r|Srr,. 



(23) 



where the second term on the right side of the equation 
is the minimum allowed value of /i, resulting from the 
finite size of the star. Note that in Papers I and II, we 
neglected to include this term in calculating /i, which 
resulted in more marked shadowing in the disk. 

For the perturbed disk, the points s^^j, s^+ij, 
s^+i^j+i, and Sij-^i are not necessarily co-planar. 
In this case, we subdivide the surface element into 
four triangles defined by the midpoint and two adja- 
cent points, i.e. (s^j, Si+i,^, s^n), (si+i,^, s^+i^^+i, s^n), 
(si+ij+i,Sij+i,Sm), and (sij+is^^^, s^n), with normal 
vectors rii, n2, 113, and 114, respectively. The angles of 
incidence are calculated as /i^ = Sm • + 4i?*/(37r|sm|), 
and the contributions to radiative heating are summed 
over each triangle individually, with Fi^r calculated at Sm 
for all of them. 

In order to avoid a discontinuity of temperature above 
the surface of the disk, the temperature above the sur- 
face is taken to be 5(rs,/i) where /i is now the angle of 
incidence to the surface of constant r.: 



IrllVr, 



(24) 



Surface elements outside the simulation box are as- 
sumed to follow the same structure as the initial con- 
ditions, i.e. unperturbed by a planet. Regions with 
101 ^ 0max are approximated as "strips" extending from 
(/>max to 7r/2 and — (/)max to — 7r/2 . The surfaces inte- 
rior to rmin and exterior to rmax are approximated as 
power laws {h ex r^) and integrated over as a series of 
infinite-length strips extending in the <p direction. 

3.3. Initial Conditions 

Since we will be calculating the perturbed disk itera- 
tively to take into account feedback between temperature 
and density perturbations, the initial disk itself must be 
stable to iterative feedback. 

We calculate the initial conditions for a slice in the 
middle of a box twice the size in (/)-space. That is, the 
slice is at = in a box ranging from — ^max to ^max 
with twice the number of grid points in and the same 
range of r and z. Outside the range of 0, the surface 
of the disk is approximated as strips from 0max to 7r/2. 
Outside the range of r, the surface is approximated as 
strips with infinite length in the (j) direction, with the 
height z varying as a power law in r, fitted to the 8 
points at the ends of the box. 

The unperturbed (i.e. without planet) structure of the 
disk is calculated iteratively as described above, but since 
azimuthal symmetry is assumed, this computation is sig- 
nificantly shorter than the computation of the perturbed 
disk structure. Thus, we iterate until the rms change in 
disk height is less than 10~^, typically less than about 
200 iterations. 

We examine planets at a = 1, 2, 4, and 8 AU with 
masses of 10, 20, and 50 Mq. The sizes and positions 
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density (g cm ^) temperature (K) 




Fig. 2. — Initial density (left) and temperature (right) profiles for the planet masses and distances studied in this paper. The lower set of 
plots are blow-ups of the area delimited by grey lines. Positions of planets at 1 (red), 2 (magenta), 4 (blue), and 8 (green) AU are indicated 
by arrows along the r-axis. The density and temperature contours are indicated by dotted, dashed, and solid curves for 10 M0, 20 Mq, 
and 50 M0 planets, respectively. 



of the simulation boxes for each of these parameters is 
indicated in Figure [2l The hmits of the simulation boxes 
in the r direction are indicated by straight vertical lines: 
dotted, dashed, and solid for 10 Mq, 20 Mq, and 50 
M0 planets, respectively. The boxes span from rmin = 
a — 12.5rHiii to rmax = a -\- 12.5rHiii in r, and from to 
^max ~ where Zs is the height of the unperturbed 
disk surface at r = a. 

The left set of plots in Figure [2] show the density con- 
tours of the initial conditions for each set of planet pa- 
rameters, dotted, dashed, and solid curves for 10 Mq, 20 
M0, and 50 planets, respectively. The lower plot is a 
blow-up of the grey-outlined box in the upper plot. The 
density contours show a good amount of overlap over the 
range of parameter space. 

The right set of plots in Figure [2] is similar to the 
left set, except it shows temperature rather than den- 
sity. The temperature contours do not overlap as well 
as the density contours, particularly at larger radii. At 
the planets positions, the temperature contour converge 
fairly well, but tend to diverge at the boundaries of the 
simulations. This results from the sensitivity of the tem- 
perature structure to the density structure of the disk. 
At small distances, the temperatures are domniated by 
viscous heating, so the shape of the surface is less rele- 
vant. However, viscous heating drops off more rapidly 
than stellar irradiation heating, so at large distances, 
slight deviations in the calculation of the surface can lead 
to changes of several degrees in the vertical temperature 
structure of the disk. The discrepancies in temperature 



may be overcome by increasing the size of the simulation 
box. However, this leads to a corresponding increase in 
the computation time. Since the difference in tempera- 
tures tends to be on the order of just a few degrees and 
since we are interested in the changes in disk structure 
due to the influence of a planet in the disk, at this time 
we will adopt the initial conditions as presented here. 

4. RESULTS 

The planet is instantaneously inserted into the initial 
conditions and the resulting density and temperature 
perturbations are calculated. The density and tempera- 
ture are then iteratively recalculated under the assump- 
tion of hydrostatic equilibrium for 10 iterations. Because 
the planet is instantaneously inserted into the disk and 
the assumption of hydrostatic equilibrium, there is no 
time scale associated with each iteration. Rather, the 
goal is to iterate until steady state is achieved. When 
evolution is discussed in this paper in the context of 
our disk-planet models, we refer to the changes in disk 
structure over successive iterations rather than a time 
sequence. 

4.1. Self- Consistency: Effect of Iteration 

The iterative calculation was implemented to achieve 
self-consistency between the temperature and density 
profiles of the disk under hydrostatic equilibrium. This 
self-consistency proves to be important to determining 
the structure of the disk. As an example, we examine 
the case of a 20 planet at 4 AU. In Figure [3] we 
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Fig. 3. — Evolution of the surface (left) and photosphere temperature (right) for a 20 earth mass planet at 4 AU. The iteration number is 
indicated at the lower right in each individual plot, (left) The deviation of the surface from the initial conditions is displayed as contours, 
in units of AU. The grey circle is a Hill radius in size, centered at the planet's location, (right) Contour plots of temperatures in the 
photosphere, in kelvins. 

dramatic by comparison, on the order of 10% or more. 

With successive iterations, the perturbation to the sur- 
face grows in area and deepens. This is seen in both the 
surface contours and photospheric temperatures. This 
is because the region in shadow cools and compresses, 
deepening the shadow further stih. Beyond the shadow, 
where the disk rises above the shadow, material heats 
and expands, causing this material to rise still further. 
The growth of the perturbation is limited both because 
of the differential rotation of disk material and because 
of the intrinsic radial temperature variation of the disk. 
The differential rotation of the disk means that mate- 
rial more radially distant from the planet moves past the 
planet faster, so it has less time to heat or cool as it 
passes by the planet. The intrinsic temperature struc- 



show the evolution of the surface perturbation (left) and 
corresponding evolution of the temperature in the pho- 
tosphere of the disk (right). The sequence of plots goes 
from left to right, top to bottom. The grey shaded circle 
shows the projected size of the Hill sphere, looked down 
on the disk along the z-axis. 

Contours in the left panel show the fractional deviation 
of the surface at the given iteration number from the ini- 
tial conditions, spaced at intervals of 0.002. Contours in 
the right panel show temperatures in the photosphere of 
the disk, where the photosphere is defined to be where 
the optical depth using the Rosseland mean opacity inte- 
grated along the 2:-axis from 2: ^ oo toward z = equals 
2/3. Although the perturbation is quite subtle, less than 
2%, the change in the temperature at the photosphere is 
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3.0 3.5 4.0 4.5 5.0 3.0 3.5 4.0 4.5 5.0 
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Fig. 4. — Evolution of temperatures in a cross-section of the disk 
at (/) = for a 20 earth mass planet at 4 AU. The Hill radius and 
position of the planet are indicated by the grey half-ellipse. The 
temperatures are indicated by contours, with the iteration number 
indicated at the lower right in each individual plot. 

ture of the disk means that inwards (outwards) of the 
planet, the disk is hotter (cooler), limiting the growth 
of the shadow (brightened region) in that direction. It- 
erations 8 — 10 show very little difference between each 
other, indicating that the density and temperatures have 
reached self-consistency. 

Figure [H shows another view of the effect of iteration 
on the temperature structure of the disk. Here we plot 
the temperature cross-section of the disk at ^ = 0, the 
location of the planet. The size of the Hill sphere is rep- 
resented by the grey ellipse (the r and z-axes are not 
scaled to each other). The surface layers are relatively 
unaffected by the the changing density structure, largely 
because this region is optically thin, and the dependence 
of temperature on angle of incidence is small. The mid- 
plane layer is relatively unaffected as well. This is due 
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Fig. 5. — Evolution of densities in a cross-section of the disk 
at = for a 20 earth mass planet at 4 AU. The Hill radius 
and position of the planet are indicated by the grey half-ellipse. 
The density varies subtly with successive iterations, so only density 
contours after the first iteration (dotted lines) and tenth iteration 
(solid lines) are shown. 

to several effects: viscous heating is greatest at the mid- 
plane; the angular sizes of the shadowed/brightened re- 
gions are small; and the disk is optically thick. The in- 
termediate layers, change in temperature structure quite 
a bit. The development of a cooled region directly above 
the planet and the heated region outward from the planet 
are quite clearly evident. 

Meanwhile, the density structure changes only subtly. 
In Figure [SI we plot the density contours for the ^ = 
slice of the simulation box after 1 (dotted lines) and 10 
(solid lines) iterations. Over the 10 iterations, the den- 
sity contours do change, but not as dramatically as the 
temperature perturbations. This emphasizes the impor- 
tance of the detailed density structure to the calculation 
of radiative transfer in the disk. 

4.2. Variations with Planet Mass and Distance 

We investigate the changes to the disk temperature 
structure as planet mass and distance vary. We exam- 
ine planets with masses of 10, 20 and 50 at 1, 2, 
4, and 8 AU. Figure [6] summarizes the changes in pho- 
tosphere tempeartures for this suite of planet parame- 
ters. The contours in each of these plots indicate the 
absolute temperatures, while the greyscale reflects the 
fractional temperature variation from the initial condi- 
tions, dark for cooling and lighter for heating. In gen- 
eral, the cooled region is centered at the planet's po- 
sition, while the headed region is outside the planet's 
orbit. The asymmetry in the perturbations are a result 
of including the differential rotation of the disk and a 
finite time scale for cooling/heating as a parcel of gas 
passes through a shadow/brightening. The sizes of the 
perturbations roughly scale both with planet mass and 
distance. 

In Figs.FZl fTOl we show radial cross-sections of the tem- 
perature structure at (from left to right) 1 Hill radius 
downstream from the planet, at the planet position, and 
1 Hill radius upstream of the planet, for the planet pa- 
rameters examined in this study. The solid lines show 
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Fig. 6. — Temperatures in the photospheres of disk models with planets after 10 iterations. From top to bottom: planets at 1 AU, 2 AU, 
4 AU, and 8 AU. From left to right: planets with mass 10 Mq, 20 Mq, and 50 Mq. The region above the Hill sphere is blacked out. Solid 
contours indicate temperatures in kelvins. The shading shows the fractional deviation in temperature from the initial disk model, so that 
black is AT/T = -0.2 and white is AT/T = +0.2 



the temperature contours, while the grey scale shows the 
aniount of deviation from the initial conditions. Figs. [7j 
[HI [HI and [To] show planets at 1, 2, 4, and 8 AU, respec- 
tively. For reference, the locations of the photospheres 
are indicated as black dotted lines. Also shown are the 
surface of the disk (white dot-dashed line) and thermal 
scale height (black dashed line). Both the photosphere 
and surface are multiple scale heights above the disk, 
except just above the Hill sphere. Temperature pertur- 
bations tend to be greatest in the upper layers of the 
disk, between the photosphere and the surface. Again, 



this is because of viscous heating at the midplane as well 
as optical depth effects. 



5. DISCUSSION 



5.1. 



Comparison to Previous Models 

We now compare our results to those of Papers I and II. 
The results presented there were for planets around a star 
of 0.5 M0, while those presented here are for a 1 Mq star. 
In order to make a fair comparison, we have recomputed 
the model presented in Papers I and II for the stellar 
parameters used in this paper. We will henceforth refer 
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Fig. 7. — Temperature cross-sections of the disk in the vicinity of a planet at 1 AU from a 1 M© star at the indicated azimuthal angles, 
corresponding to (from left to right) 1 Hill radius downstream of the planet, at the planet position, and 1 Hill radius upstream of the 
planet. The location and size of the Hill sphere are indicated by the white hashed area. From top to bottom, planet masses are 10 Mq, 20 
M0 and 50 Mq. The contours show absolute temperatures while the color scale show the fractional temperature difference from the initial 
disk model, with a range of AT/T G [-0.2, +0.2]. 



to this set of models as the J-CS model. 

In Figure [m we compare the results from this paper to 
J-CS for the case of a 20 planet at 4 AU from its host 
star. The top plot shows the radial temperature cross 
section through the planet's position for J-CS and the 
bottom shows the same cross-section for the new model. 
In the J-CS model, the disk is treated as locally plane- 
parallel, whereas in the new model it is not. This is 
reflected in the shape of the temperature structures of 
the two models. Another difference is that the surface 
was calculated as the contour of constant density in J-CS. 



This contour is shown as a red dot-dashed line. When 
this contour is shadowed, the surface is instead taken to 
be the line of sight to the star, indicated by the solid red 
line. In contrast, the new model calculates the surface to 
be where the optical depth to the star using the Planck 
mean extinction, Xp, is 2/3. 

In J-CS, the temperature perturbations are more dra- 
matic than in the new model. This can also be seen in 
Figure [121 which shows the temperatures in the photo- 
sphere compared side- by-side. The plots show tempera- 
ture contours and fractional temperature variations, with 
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Fig. 8. — Same as Figure [T] for planets at 2 AU. 



the older model represented on the left. The right plot 
in this figure is the same as the one in Figure [6l The 
spatial sizes and color scalings are the same between the 
two plots. It is evident from both Figs. [11] and [12] that 
the temperature perturbations are much greater in J- 
CS. The reason for this is primarily because of the way 
in which the surface was calculated. As noted in Paper 
I, the disk temperature structure is quite sensitive to /i, 
the cosine of the angle of incidence of stellar irradition at 
the surface. In other words, the slope of the surface de- 
termines the amount of heating from stellar irradiation. 
Small deviations at the surface can lead to large temper- 
ature perturbations at the photosphere of the disk. 
In the JC-S model, /i ^ in the shadowed region. 



which means that zero fiux was propagated from shad- 
owed surface elements to disk material below them. 
When the density contour rose back above the shadow, 
the transition from shadow to the illuminated region 
was quite sharp, with relatively large values of fi yield- 
ing large amounts of heating. In the new model, since 
the surface is calculated by integrating the optical depth 
along lines of sight to the star, the transitions are 
smoother and values of fi do not vary as greatly. Also, 
as shown in equation ([23]) , we set a minimum value for fi 
based on the finite size of the star, so shadowing can not 
result in as significant cooling as seen in J-CS. Yet an- 
other effect is the radial variation in stellar fiux, which 
was not taken into account in the plane-parallel J-CS 
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Fig. 9. — Same as Figure [T] for planets at 4 AU. 



model. The disk is hotter/cooler closer to/farther from 
the star, which offsets some of the effects of shadow- 
ing/illumination. Although the J-CS model was not it- 
erative while the new model is, the changes in the calcu- 
lation of the surface and /i are sufficient to diminish the 
temperature perturbations near the planet. 

Another difference between J-CS and this work is the 
locations of the heated and cooled regions. In J-CS, these 
regions are nearly symmetric with respect to the planet 
position. The sizes and shapes of the heated and cooled 
regions are similar, with the cooled region inward of the 
planet and the heated region outward of it. In the cur- 
rent work. The cooled region in centered on the planet's 
position with the heated region outward of the planet. 



5.2. Magnitude of Temperature Variations 

In Paper II, we found a correlation between the maxi- 
mum and minimum photosphere temperature variations 
and the Hill radius of the planet. We now perform the 
same analysis on the new model. 

In Figure [131 we plot the maximum and minimum frac- 
tional temperature variations in the photospheres versus 
the ratio of the Hill radius to thermal scale height of the 
disk {ruiw/h). Points above/below AT/T = represent 
maxima/minima. The symbols denote at what distance 
the planet is: triangles for 1 AU, diamonds for 2 AU, 
squares for 4 AU, and asterisks for 8 AU. The results for 
J-CS are plotted connected by dashed lines. As noted in 
Paper II, the maxima and minima lie nearly on the same 
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Fig. 10. — Same as Figure [3 for planets at 8 AU. 

curves. The upturn in minimum temperatures at smaller 
distances is due to the rise in viscous heating closer to 
the star, which sets the absolute temperature minimum. 

The points connected by dotted lines in Figure [13] are 
temperature maxima/minima in the photospheres after 
one iteration, i.e. without recalculating the density given 
the new temperature structure. These points similarly 
lie on nearly the same curves, except showing a much 
smaller temperature perturbation than JC-S. The rea- 
sons for the smaller temperature variation was detailed 
in go 

The effect of including iterations can be seen by com- 
parison to the points connected by solid lines in Figure 
[T3l In all cases, the fractional temperature variations 



increase. However, the amount of change is not consis- 
tent with planet mass or distance. At 1 AU (triangles), 
the change in temperature minima is nearly negligible. 
This can be explained by viscous heating, which halts 
the growth of the cooled, shadowed region. The temper- 
ature maxima do increase, but not as much as at larger 
distances. This is because the growth of the brightened 
region is coupled to the growth of the shadowed region: 
as the shadowed region grows, more disk material behind 
the shadow is exposed and illuminated. Since shadow 
growth is limited at 1 AU, so is the growth of the bright- 
ened region. Going from 2 to 4 to 8 AU, viscous heating 
ceases to be important. The temperature maxima all lie 
close to the same curve at those distances, and the tem- 
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Fig. 11. — The thermal structure of the disk in the vicinity of a 
20 M0 planet at 4 AU from a 1 star, as calculated for (top) 
J-CS and (bottom) this paper. The temperatures are scaled to 
the same colors, and the r and z axis have the same scaling. The 
hashed semi-circle shows the size and location of the Hill sphere. 
The photosphere is indicated by white dotted lines and the thermal 
scale height is indicated by black dashed lines. In the top plot, 
the red dot-dashed line shows the isodensity contour equal to the 
density at the disk surface. Regions below the solid red line in the 
top plot are considered to be in shadow. In the bottom plot, the 
red dot-dashed line is the surface as calculated from line-of-sight 
optical depth to the star. 

perature minima at 4 and 8 AU lie nearly on the same 
curve in Figure [T3l 

For larger planets, those whose Hill radii approach the 
disk thermal scale height, the amount of heating in the 
current model approaches that seen in J-CS. However, 
the temperature minima never drop that low. This can 
be attributed to the minimum value of ji that has been 
adopted for the new model. There is such upper limit 
that has been imposed however, hence the greater heat- 
ing. 



the photosphere that are heated or cooled above or be- 
low a certain threshold. We call these regions hot or cold 
"spots." Note that this is different from the definition of 
a "spot" in Paper H, which considered regions heated 
or below 170 K, in the context of ice formation and the 
snow line. Here, we consider the deviation in tempera- 
ture from the unperturbed value, rather than above or 
below a fixed threshold. A separate paper addressing the 
effect of these tem perature pert urbations on the snow 
line is addressed bv I Jang-Condell et al.l (2007). 

We define the area of a spot to be the area of the 
photosphere that is heated or cooled to at least half the 
maximum temperature deviation, excluding the area just 
above the Hill sphere. These spot areas are plotted as 
filled triangles (hot) or squares (cold) in Figure [TH con- 
nected by dashed lines, black for the current models, 
grey for J-CS. The temperature deviations for the 10 
M0 planet at 8 AU are so small that the area of a cold 
spot is not well-defined. The spot masses are generally 
larger for our results than for the JC-S models. The sizes 
of cold spots tend to grow faster with increasing planet 
mass than the sizes of hot spots. 

Let us also consider temperature perturbations above 
or below a threshold of 10%. The area of these spots 
versus the area of a circle with radius rniii is plotted as 
open triangles (hot) or squares (cold) in Figure [TH con- 
nected by dotted lines, black for the current models, grey 
for J-CS. Lines connected to points below the edge of the 
graph indicate lack of a spot. Note that at 1 AU, there 
are no hot spots above a 10% temperature deviation, con- 
sistent with Figure [T3l The J-CS models produce hot at 
cold spots for planet parameters. At all distances, 10 M0 
planets are too small to create either hot or cold spots 
at the 10% level in the new models. Where cold spots 
do exist, they are generally larger than those produced 
in J-CS. However, the hot spots are generally smaller. 

While the temperature deviations predicted from the 
current models are not as great as those predicted in 



5.3. Sizes of Thermally Perturbed Regions 

Although the variation in temperature appears to be 
less in the present model than that seen in J-CS, Figs.fTTl 
and [12] suggest that the spatial scale of the perturbations 
might be larger. To check this, we calculate the areas of 
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Fig. 12. — Temperature perturbations in the disks photosphere 
above an embedded planet, (left) as calculated in the J-CS model 
and (right) using the methods presented in this paper. Figures 
are scaled to the same spatial size, and colors represent fractional 
deviation from the unperturbed disk temperature. Contours show 
temperatures. Since the J-CS calculation assumed a plane-parallel 
disk, the temperature contours trace the color variations. In the 
calculation presented in this paper, the temperatures are domi- 
nated by the radial variation of the unperturbed disk. 
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Fig. 13. — Maximum and minimum fractional temperature vari- 
ation in the photosphere versus ^Hill/^- Symbols indicate planet 
distance as follows: triangles at 1 AU, diamonds at 2 AU, squares 
at 4 AU, and asterisks at 8 AU. Dashed lines show the J-CS model. 
Dotted lines show the current model, after only one iteration. Solid 
lines show the current model, after 10 iterations. 
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Fig. 14. — Areas of heated/cooled regions in the photosphere 
in current model (black) compared to results of previous model 
(grey). Horizontal axis is the area of a circle with radius of rj^in. 
Planets with 10, 20 and 50 Me at 1, 2, 4, and 8 AU are shown. The 
distances for each set of planets are indicated by the long-dashed 
boxes. Heated and cooled regions are represented by triangles and 
squares, respectively. Filled symbols connected by dashed lines 
show the areas heated/cooled beyond half the maximum temper- 
ature deviation in the photosphere. Open symbols connected by 
dotted lines show the areas that are heated/cooled beyond 10% 
above or below the unperturbed temperature. The diagonal solid 
lines show a linear slope. 

J-CS, they are larger in spatial extent. This may have 
consequences for observability of this phenomenon. This 
topic will be addressed in a companion paper. 

6. CONCLUSION 

This paper presents a model for calculating the density 
and temperature perturbations imposed on a protoplan- 
etary disk by an embedded protoplanet. The basic radia- 
tive transfer model is adopted from Papers I and II, but 
a number of improvements have been made on that work 
such as density-temperature self-consistency and elimi- 
nating the assumption of a plane-parallel system. 

We have shown that self-consistently calculating the 
temperature and density significantly increases the ef- 
fect of planet perturbations by means of postive feed- 
back, where shadowed regions cool and contract and 



brightened regions heat and expand. This demonstrates 
the importance of self-consistency when calculating disk 
structure with radiative transfer. While it has al- 
ready been acknowledged that stellar irradiation heat- 
ing is i mportant in setting overall disk structur e (e.g . 
Chiang Goldreich 19 97; D ^Alessio et al. 1999^ fT998l : 



Dullemond fc Dominikll2QQ4l ). this work shows that it is 



important for considering local perturbations on the disk. 

Another important result of this paper is that the tem- 
perature structure of the disk is extremely sensitive to the 
angle of incidence of stellar irradiation at the surface. 
The precise determination of the surface of the disk is 
critical to accurately calculating the temperature struc- 
ture of the disk. In previous work we had assumed that 
the surface of constant density was a sufficiently good ap- 
proximatation for the surface, but in this work we show 
that that overestimated the temperature perturbation to 
the disk. 

We have omitted some important physics in this calcu- 
lation. We do not include heating from accretion onto the 
planet. We consider the embedded planet to act simply 
as a gravitational point mass, and focus only on the ver- 
tical component of gravity. We do not account for non- 
linearities such as spiral density waves. All these effects 
are more adequately addressed using a three-dimensional 
hydrodynamic simulation of a planet embedded in a disk 
of gas. However, since simulations of this sort focus on 
the bulk flow of gas at the midplane, they typically have 
insufficient resolution above the midpl ane to accurately 
calcul ate the s urface of the dis k (e.g; . [Bate et al."2 QQ3|: 
Papaloizou et al. 2004; Klahr KlevI P006; Oishi et^ 
I2OO7I ). Calculation of radiative transfer, even without 
iterating for self-consistency, is very computationally in- 
tensive. To do this iteratively and coupled with three- 
dimensional hydrodynamics is challenging, but is the 
next logical step to improving the accuracy of our re- 
sults. 
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